% By Martin M. Andreasen

function [model,errorMes] = perturbationDSGEmodel(params,appOrder,setupModel,maturities,MexOn)

% Computing the steady state and checking that nf = 0
[auxOut,nf,errorMes] = numF(params);
params = auxOut.params;
if any(abs(nf)> 1D-8) 
    disp('Displaying nf:')
    disp(nf)
end
if errorMes ~= 0
    model = NaN;
    errorMes = 1;
    return;
end

% The dimension of the model
eta = auxOut.eta;
ne  = auxOut.ne;
nx  = auxOut.nx;
ny  = auxOut.ny;

% Ensure same ordering of params as in symparams
paramsOrdered = struct();
for i=1:length(setupModel.symparams)
    paramsOrdered.(char(setupModel.symparams(1,i))) = params.(char(setupModel.symparams(1,i)));
end
     
% Choose solution algorithm: 
% algo='vectorize'; % standard vectorization - good only for small models
%algo='dlyap'; % Hessenberg-Schur algorithm - good for large models, but
%                 requires to have the control system toolbox
algo='gensylv'; % Kamenik (2005) algorithm - recommended. If the MEX file gensylv.mexw64 does not work on your system, download dynare from www.dynare.org, search for gensylv and add the MEX file to the Perturbation folder.

% Solve the perturbation loadings
derivs = solve_dsge(setupModel,struct2array(paramsOrdered),...
    auxOut.momEps,eta,auxOut.xssTrans',auxOut.yssTrans',appOrder,algo,MexOn);

if isempty(derivs)
    model = NaN;
    errorMes = 1;
    return;
end

%% Saving the output in model
model.labelx           = auxOut.labelx;
model.labely           = auxOut.labely;
model.g0               = auxOut.yssTrans;
model.h0               = auxOut.xssTrans;
model.derivs           = derivs;
model.gx               = derivs.gx(1:ny,1:nx);
model.hx               = derivs.hx(1:nx,1:nx);
if appOrder > 1
    tmp_gxx            = reshape(derivs.gxx,ny,nx+1,nx+1);
    model.gxx          = tmp_gxx(1:ny,1:nx,1:nx);
    model.Gxx          = 1/2*reshape(model.gxx,ny,nx^2);
    model.gss          = tmp_gxx(1:ny,nx+1,nx+1);
    tmp_hxx            = reshape(derivs.hxx,nx,nx+1,nx+1);
    model.hxx          = tmp_hxx(1:nx,1:nx,1:nx);
    model.Hxx          = 1/2*reshape(model.hxx,nx,nx^2);
    model.hss          = tmp_hxx(1:nx,nx+1,nx+1);
else
    model.gxx          = zeros(ny,nx,nx);
    model.Gxx          = zeros(ny,nx^2);
    model.gss          = zeros(ny,1);
    model.hxx          = zeros(nx,nx,nx);
    model.Hxx          = zeros(nx,nx^2);
    model.hss          = zeros(nx,1);
end

if appOrder > 2
    tmp_gxxx           = reshape(derivs.gxxx,ny,nx+1,nx+1,nx+1);
    model.gxxx         = tmp_gxxx(1:ny,1:nx,1:nx,1:nx);
    model.Gxxx         = 1/6*reshape(model.gxxx,ny,nx^3);
    model.gssx         = squeeze(tmp_gxxx(1:ny,nx+1,nx+1,1:nx));
    model.gsss         = squeeze(tmp_gxxx(1:ny,nx+1,nx+1,nx+1));
    tmp_hxxx           = reshape(derivs.hxxx,nx,nx+1,nx+1,nx+1);
    model.hxxx         = tmp_hxxx(1:nx,1:nx,1:nx,1:nx);
    model.Hxxx         = 1/6*reshape(model.hxxx,nx,nx^3);
    model.hssx         = squeeze(tmp_hxxx(1:nx,nx+1,nx+1,1:nx));
    model.hsss         = squeeze(tmp_hxxx(1:nx,nx+1,nx+1,nx+1));
else
    model.gxxx         = zeros(ny,nx,nx,nx);
    model.Gxxx         = zeros(ny,nx^3);
    model.gssx         = zeros(ny,nx);
    model.hxxx         = zeros(nx,nx,nx,nx);
    model.Hxxx         = zeros(nx,nx^3);
    model.hssx         = zeros(nx,nx);
end

%% Step 2: Solving for bond prices by POP
% Solve for bond prices
[nM,nMxp,nMyp,nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp] = numDerivM(params);
Max_maturity = max(maturities);
vectorMom3   = zeros(1,ne);
ny_1bond     = find(strcmp(auxOut.labely,'$p1_t$'));
solveRiskTermsSmallModel = 1;
if MexOn == 1
    [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3_LogMex(nM,nMxp,nMyp,...
    nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
    model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
    model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
    eta,vectorMom3,int32(Max_maturity),int32(ny_1bond),int32(appOrder),int32(solveRiskTermsSmallModel),int32(ny),int32(nx),int32(ne));
else
    [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3_Log(nM,nMxp,nMyp,...
        nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
        model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
        eta,vectorMom3,Max_maturity,ny_1bond,appOrder,solveRiskTermsSmallModel);
end

%% Compute conditional expectations of the short rate
y_select = find(strcmp(auxOut.labely,'$r_t$')); % The short rate is stored at position at element 2
% The arrays contain expectations from 1,2,..., Max_maturity into the future
maxMatShortRate = 4;
[rx,rxx,rss,rxxx,rssx,rsss] = CondMoments_Mom3_log(model.gx,model.gxx,model.gss,...
    model.gxxx,model.gssx,model.gsss,...
    model.hx,model.hxx,model.hss,...
    model.hxxx,model.hssx,model.hsss,eta,...
    vectorMom3,y_select,maxMatShortRate,appOrder);

% Adding bond yields to model
model.by0    = -bsxfun(@ldivide,(1:Max_maturity)',log(P(:,:))); 
model.byx    = -bsxfun(@ldivide,(1:Max_maturity)',Px(:,:));
model.byxx   = -bsxfun(@ldivide,(1:Max_maturity)',Pxx(:,:,:));
model.Byxx   = -0.5*reshape(bsxfun(@ldivide,(1:Max_maturity)',Pxx(:,:,:)),Max_maturity,nx^2);
model.byss   = -bsxfun(@ldivide,(1:Max_maturity)',Pss(:));
model.byxxx  = -bsxfun(@ldivide,(1:Max_maturity)',Pxxx(:,:,:,:));
model.Byxxx  = -1/6*reshape(bsxfun(@ldivide,(1:Max_maturity)',Pxxx(:,:,:,:)),Max_maturity,nx^3);
model.byssx  = -bsxfun(@ldivide,(1:Max_maturity)',Pssx(:,:));
model.bysss  = -bsxfun(@ldivide,(1:Max_maturity)',Psss(:));
model.nYields= Max_maturity;
model.maturities = maturities;
% Adding conditional expectations to model
model.matConShortRate = 1:maxMatShortRate;
model.r0     = repmat(model.by0(1,1),maxMatShortRate,1);
model.rx     = rx;
model.rxx    = rxx;
model.Rxx    = 1/2*reshape(rxx,maxMatShortRate,nx^2);
model.rss    = rss;
model.rxxx   = rxxx;
model.Rxxx   = 1/6*reshape(rxxx,maxMatShortRate,nx^3);
model.rssx   = rssx;
model.rsss   = rsss;
model.appOrder= appOrder;
model.ny     = ny;
model.nx     = nx;
model.ne     = ne;
model.nx1    = nx-ne;
model.eta    = eta;
model.params = params;
